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1  •  INTRODUCTION 

The  computation  of  a  detailed  geoid,  or  of  a  detailed  gravity  poten¬ 
tial  field,  in  limited  areas,  especially  in  mountainous  regions,  has  not 
been  very  much  in  the  focus  of  attention  recently.  There  may  be  various 
reasons  for  this. 

For  two  decades  now,  gl obal  geoid  determinations,  either  from  satel¬ 
lite  data  or  from  a  combination  of  satellite  and  gravimetric  data  have  been 
in  the  center  of  interest  (cf.  Lerch,  et  al.,  1979;  Reigber,  et  al . ,  1983; 
Rapp,  1981).  Even  (almost)  purely  gravimetric  global  geoids  have  been 
successfully  computed  (cf.  March  and  Chang,  1979). 

Over  the  oceans,  the  geoid  is  now  known  to  an  accuracy  of  perhaps  +_1 
or  +2m,  due  to  satellite  altimetry.  Unfortunately,  satellite  altimetry 
does  not  work  over  land  areas.  The  classical  method  for  a  detailed  geoid 
determination  on  the  continents  is  the  gravimetric  method,  in  spite  of  the 
fact  that  it  is  severely  handicapped  by  lack  of  an  adequate  gravity 
coverage  (or  lack  of  information  on  such  a  coverage).  Thus  we  have  the 
paradoxical  situation  that  on  the  oceans,  long  a  stepchild  of  geodesy,  the 
geoid  is  now  in  general  known  much  better  than  on  the  continents. 

Still,  the  gravimetric  method  has  continued  to  fascinate  theoreticians 
because  it  gives  rise  to  very  interesting  and  deep  mathematical  problems, 
related  to  the  geodetic  boundary-value  problem,  or  problem  of  Molodensky 
(cf.  Moritz,  1980,  Part  D).  In  the  recent  years,  the  combination  of  satel¬ 
lite  altimetry  on  sea  and  gravimetry  on  land  has  led  to  another  interesting 
boundary-value  problem,  the  altimetry-gravimetry  boundary-value  problem 
(cf.  Sanso,  1983). 

These  enormous  practical  and  theoretical  developments  concerning 
global  satellite  and  gravimetric  gravity  field  determination  have  somewhat 
overshadowed  the  determination  of  detailed  geoids  in  smaller  areas. 
Especially  in  mountainous  regions,  such  local  geoid  determinations  are 
difficult.  The  gravimetric  method  does  not  work  very  well  in  high 
mountains.  The  astrogeodeti c  method,  using  astronomical  observations  of 
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latitude  and  longitude,  does  work  well  there,  but  is  considered  time- 
consuming  and  somewhat  old-fashioned,  perhaps  also  because  working  during 
the  night  is  not  very  popular  nowadays.  An  appropriate  use  of  gravity  and 
astrogeodetic  data  in  high  mountains  must  involve  some  topographic- 
isostatic  reduction,  which  is  also  sometimes  considered  old-fashioned. 
Furthermore,  the  theory  behind  the  astrogeodetic  method  is  not  nearly  as 
attractively  difficult  as  the  theory  of  Molodensky's  problem.  Last  but  not 
least,  high-mountain  areas  are  exceptional  and,  apart  from  such  countries 
as  Switzerland  and  Austria,  are  frequently  regions  of  little  economic 
interest.  For  these  and  similar  reasons,  the  main  stream  of  geodetic 
practice  and  theory  has  flown  with  grand  indifference  around  high 
mountains,  ignoring  such  trivial  obstacles. 

Still,  a  country  such  as  Switzerland  has  made  a  virtue  out  of  neces¬ 
sity  and  has  traditionally  been  very  active  in  local  astrogeodetic  geoid 
determination  (Elmiger,  1969;  Gurtner,  197ft;  Gurtner  and  Elmiger,  1983). 
Recently,  Austria  has  followed  up  (Lichtenegger,  et  al.,  1983).  It  has 
been  found  that,  even  besides  the  problem  of  getting  the  required 
observations,  the  underlying  theory  is  not  so  trivial  as  one  might  think 
ind  shows  quite  interesting  features. 

Concerning  measurements,  astronomical  observations  have  again  proved 
very  feasible  in  mountains;  see  the  articles  by  Erker,  Bretterbauer , 
Lichtenegger  and  Chesi  in  Chapter  ?.  of  (Lichtenegger,  et  al.,  1983).  The 
main  advantages  of  astrogeodetic  versus  gravimetric  data  for  local  geoid 
determination  in  mountain  regions  may  be  summarized  as  follows: 

1.  It  is  sufficient  to  have  astrogeodetic  deflections  of  the  verti¬ 
cal  in  the  region  of  geoid  determination;  no  data  are  needed  outside  that 
region  as  they  would  be  in  the  gravimetric  method. 

7.  Errors  in  the  topographic  height  have  significantly  less  influ¬ 
ence  on  deflections  than  on  gravity  data.  Thus  a  relatively  crude  terrain 
model  will  be  sufficient  for  the  use  of  astrogeodetic  data. 

As  a  matter  of  fact,  the  two  types  of  data  are  not  mutually  exclusive; 
an  optimal  geoid  determination  will  combine  astrogeodetic  deflections  of 
the  vortical,  gravity  anomalies,  and  possibly  data  of  other  type.  A  suita¬ 
ble  technique  for  this  purpose  is  least-squares  collocation. 
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From  the  observational  point  of  view  it  is  interesting  to  note  that 
inertial  surveying  techniques  will  be  able  to  furnish  deflections  of  the 
vertical  and  gravity  anomalies  rapidly  and  with  sufficient  accuracy  for 
many  purposes. 

Let  us  finally  try  to  give  a  list  of  various  methods  of  geoid  de¬ 
termination: 

conventional  satellite  techniques  (doppler,  laser,  etc.) 

satel 1 ite-to-satel 1 ite  tracking 

satellite  gradiometry 

satellite  altimetry 

aerial  gradiometry 

gravimetry 

astrogeodesy 

As  a  general  rule,  these  methods  are  listed  in  such  a  way  as  to  start  with 
the  most  global  and  end  up  with  the  most  local  method,  that  is,  according 
to  decreasing  global ity  or  increasing  locality.  In  general,  going  down  the 
list  also  corresponds  to  increasing  resolution  and  accuracy. 

Again  it  should  be  stressed  that  these  methods  complement  each  other 
and  should  be  combined  for  best  results. 

The  present  report  deals  primarily  with  the  lower  end  of  the  list, 
providing  a  detailed  theory  of  local  geoid  determination  in  areas  with  dif¬ 
ficult  topography.  The  role  (and  necessity)  of  topographic-isostatic  re¬ 
duction  is  investigated  in  detail.  The  computations  for  Austria  give  con¬ 
crete  numerical  results  for  questions  which  have  been  much  discussed  theo¬ 
retically,  such  as  the  difference  between  geoidal  heights  and  height  ano¬ 
malies  according  to  Molodensky  (quasigeoidal  heights),  or  the  numerical 
effect  of  analytical  continuation  from  the  earth's  surface  to  sea  level. 


Definitions 


FIGURE  1.  The  basic  geometry 


Fig.  1  illustrates  the  basic  quantities.  In  the  classical  theory,  the 
geoid  is  defined  by  its  deviation  N  from  a  reference  ellipsoid;  N  is 
the  geoidal  height.  The  geoid  is  a  level  surface  VI  =  VJ  =  const,  of  the 
gravity  potential  W  ;  the  ellipsoid  is  defined  to  be  the  level  surface  U 

=  U  =  const,  of  a  normal  gravity  potential  U  ;  the  constants  VI  and  U 

o  0  0 

are  usually  assumed  to  be  equal.  Cf.  PG,  sec.  2-13*). 

For  the  modern  theory  according  to  Molodensky  (PG,  sec.  8-3),  to  each 
point  P  of  the  earth's  surface  we  associate  a  point  Q  in  such  a  way 
that  0  lies  on  the  straight  ellipsoidal  normal  through  P  and  that 

U(Q)  =  W(P)  .  (1) 


That  is,  Q  is  defined  such  that  its  normal  potential  U  equals  the 
actual  potential  W  of  P  . 

This  corresponds  to  the  classical  relation 

u„  ,  U(Q0)  =  W(P0)  =  «0 

mentioned  above,  by  which  IJQ  is  taken  to  be  equal  to  W  ;  cf.  Fig.  1. 
By  the  same  correspondence,  the  height  anomaly  according  to  Molodensky, 


?  ■  QP 


(3) 


is  the  modern  equivalent  of  the  classical  geoidal  height, 


N  = 


Q  P 

o 


(4) 


1)  By  the  symbol  PG  we  shall  in  the  sequel  denote  the  book 
"Physical  Geodesy"  (Heiskanen  and  Moritz,  1967). 


Using  the  anomalous  potential 


I 

we  have  according  to  Bruns'  theorem 


Y  denoting  ellipsoidal  normal  gravity. 

The  points  P  form  the  geoid,  and  the  points  Q  constitute  the 
ellipsoid,  both  being  level  surfaces  (of  W  and  II  ,  respectively).  On 
the  other  hand,  the  points  P  form  the  earth's  surface,  and  the  set  of 
points  0  defines  an  auxiliary  surface,  denoted  as  tel  1 uroid  according 
j  to  R.  A.  Hirvonen.  As  a  matter  of  fact,  neither  the  earth's  surface  nor 

the  Telluroid  are  level  surfaces,  which  makes  matters  more  complicated  than 
in  the  classical  situation,  where  we  deal  with  level  surfaces. 

Following  a  suggestion  of  Molodensky,  one  could  plot  the  height 
|  anomalies  C  as  vertical  distances  from  the  reference  ellipsoid.  Thus  one 

obtains  a  geoid-like  surface,  the  quasi -geoi d ,  and  £  could  be  considered 
as  quasi -geoidal  heights.  In  contrast  to  the  geoid,  however,  the  quasi- 
geoid  is  not  a  level  surface  and  does  not  admit  of  a  natural  physical 
i  interpretation.  Therefore,  working  with  height  anomalies  c,  ,  it  is  best 

to  consistently  consider  them  quantities  referred  to  the  earth's  surface 
(vertical  distances  between  earth  surface  and  telluroid),  rather  than  using 
the  quasi -geoidal  concept. 

[  The  classical  gravity  anomaly  Ag  at  sea  level  is  defined  as 


Ag„  -  g(PJ  -  y(QJ  . 


(7) 


where  g  denotes  gravity  and  y  ,  normal  gravity.  (So  far,  g(PQ)  de¬ 
notes  the  actual  gravity  on  the  geoid;  we  are  not  yet  considering  mass¬ 
transporting  gravity  reductions.) 

Analogously  we  have  according  to  Molodensky: 

Ag  =  g(P)  -  y(Q)  • 


Generally  we  shall,  as  far  as  feasible,  use  the  subscript  "o"  to  de¬ 
signate  quantities  referred  to  sea  level,  to  distinguish  them  from 
quantities  referred  to  the  earth's  surface,  which  do  not  carry  such  a  sub¬ 
script.  For  instance,  AgQ  refers  to  sea  level  and  Ag  ,  to  the  earth's 
surface. 

Regarding  plumb-line  definition,  we  must  distinguish  three  lines  (Fig. 

1): 

(a)  The  straight  ellipsoidal  normal  QqP  , 

(b)  The  actual  plumb-1 ine  P^’P  , 

(c)  The  normal  plumb-line  P^P  . 

The  ellipsoidal  normal  is  geometrical ly  defined  as  the  straight  line 
through  P  perpendicular  to  the  ellipsoid.  The  (actual)  plumb  line  is  de¬ 
fined  by  the  condition  that,  at  each  point  of  the  line,  the  tangent 
coincides  with  the  gravity  vector  g  at  that  point;  the  plumb  line  is 
very  slightly  curved,  but  its  curvature  is  irregular,  being  determined  by 
the  irregularities  of  topographic  masses.  The  normal  plumb  line,  at  each 
of  its  points,  is  tangent  to  the  normal  gravity  vector  y  ;  it  possesses  a 
curvature  that  is  even  smaller  and  completely  regular. 

The  points  P  ,  P'  ,  and  P1'  coincide  within  a  few  decimeters, 

OO  0 

and  we  shall  not  distinguish  them  in  the  sequel.  The  reason  is  that  the 
distance,  in  seconds  of  arc,  between  Pq  and  P^ '  ,  is  much  smaller  than 
the  effect  of  plumb  line  curvature  (PG,  p.  180-181).  The  same  holds,  of 
course,  for  Q0  ,  Oq  ,  and  Qq'  . 


The  direction  of  the  gravity  vector  g  is  the  direction  of  (the 
tangent  to)  the  plumb  line.  It  is  determined  by  two  angles,  the 
astronomical  latitude  $  and  the  astronomical  longitude  A  .  Let  4  ,  A 
he  referred  to  the  earth  surface  (to  point  P  )  and  <5>o  ,  AQ  to  the 
geoid  (strictly  speaking,  to  point  P^  ).  The  differences 

64  =  4  -  4  ,  6A  s  A  -  A  (9) 

o  o 


parallel  to 
equatorial  plane 


FIGURE  2.  Curvature  of  the  plumb  line  along  a  north-south  profile 


express  the  effect  of  plumb-line  curvature  (Fig.  2).  Hence  we  have 


=  $  +  <54>  , 


A  =  A  +  cSA  . 
o 


Knowing  the  plumb-line  curvature  (  Si  ,  SA  ) ,  we  could  use  these  simple 
formulas  to  compute  the  sea-level  values  iQ  ,  Aq  from  the  observed 

surface  values  $  ,  A  . 

In  the  same  way  as  i  ,  A  are  related  to  the  actual  plumb  line,  the 
geodetic  latitude  <(>  and  the  geodetic  longitude  X  refer  to  the  straight 
ellipsoidal  normal.  The  quantities 


£  ■  «  -  i  , 


n  =  (A  -X  )cosg> 


are  the  components  of  the  deflection  of  the  vertical  in  a  north-south  and 
an  east-west  direction.  For  an  arbitrary  azimuth  a  ,  the  vertical  .de¬ 
flection  e  is  given  by 

(12) 

c  =  £cosa  +  risma  .  ' 

These  quantities  £  ,  n  ,  e  refer  to  the  earth's  surface.  Cf.  Fig.  1, 
which  shows  e 

Similarly  we  have  for  the  geoid 


=  • 


n0  =  (AQ  -  x)cos4>  , 


eQ  =  £QCOSa  +  nQsina  . 
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See  again  Fig.  1  for  eq  ,  noting  that  we  do  not  distinguish  the  normals 

in  0  and  O'1  as  we  have  mentioned  above. 

o  o 

In  addition,  we  need  the  normal  direction  of  the  plumb  line  at  the 
surface  point  P  ;  it  is  defined  as  the  tangent  to  the  normal  plumb  line 
at  P  ;  the  corresponding  latitude  and  longitude  will  be  denoted  by  $  , 

A  .  Hence  we  have 


4>  =  <j>  +  6<J>  ,  A  =  A  +  6A 


(15) 


where  6<j>  ,  SA  express  the  normal  plumb-line  curvature.  These  equations 

are  the  "normal  equivalent"  to  (10):  the  "normal  surface  values"  p  ,  X 

correspond  to  the  "actual  surface  values"  4>  ,  A  ,  and  the  ellipsoidal 

values  p  ,  A  correspond  to  the  geoidal  values  $  ,  A  .  (To  make 

0  o 

the  analogy  complete,  we  should  replace  <J>  =  ( Pq )  by  4)(P^)  »  hut  we 

have  consistently  been  neglecting  such  differences.) 

In  marked  contrast  to  the  actual  plumb-line  curvature,  it  is  very  easy 
to  compute  the  normal  curvature  of  the  plumb  line:  by  PG,  p.  196  we  have 

5cp  =  -0.17' 'hkmsin2$  ,  6A  =  0  ,  (16) 

hkr,  denoting  elevation  in  kilometers. 

Since  the  ellipsoidal  normal  and  hence  P  ,  A  are  geometrically 
defined,  we  may  call  the  quantities  (11)  "geometric  deflections  of  the 
vertical"  at  the  earth's  surface.  On  the  other  hand,  the  normal  plumb  line 
is  physically  (or  dynamically)  defined  by  means  of  the  external  gravity 
field  of  an  equi potenti al  ellipsoid.  Hence  also  P  ,  A  are  dynamically 
defined,  and  we  may  call  the  quantities  obtained  by  replacing  <t>  ,  A  by 

1  ,  X  : 


JP. 


£  =  $  -  P  > 


n  =  (  A-A  )cos0  , 


(17) 


and  (16)  we  have 


C  =  S  +  <54>  >  n  =  n 


since  6X  =  0  .  For  an  azimuth  a  we  accordingly  have 
e  =  £co sa  +  nsina 

Compare  e  and  T  in  Fig.  1,  and  note  that  in  this  figure, 
the  curvature  of  the  normal  plumb  line  for  the  azimuth  a  , 
analogous  formula 


6  =  6<|>cosa  +  (6Xcos<t>)sina  =  6<j>cosa  . 


(18) 


(19) 

6  denotes 
given  by  the 


(20) 


3. 


FIGURE  3.  Astronomical  Leveling  according  to  Helmert 
From  Fig.  3  we  take  the  well-known  differential  relation 


dN  =  -eQds  , 

denoting  the  deflection  of  the  vertical  at  the  geoid.  Integration 
between  two  points  A  and  B  yields  the  difference  between  their  geoidal 
heights: 


or  approximately, 


where  sab  denotes  the  horizontal  distance  between  A  and  B.  The  minus 
sign  is  conventional. 


FIGURE  4.  Astronomical  leveling  according  to  Molodensky 

A  corresponding  relation  to  height  anomalies  according  to  Molodensky 
is  found  as  follows  (Molodensky,  et  al.,  1962,  p.  125): 


d£ 


3S 


ds  + 


3h 


dh 


(24) 


notations  following  Fig.  4.  Since  the  earth's  surface  is  not  a  level  sur¬ 
face,  we  also  have  a  vertical  part  (  3c/3h  )  h  in  addition  to  the 

usual  horizontal  part  (  3;; /3s  )  ds  .  The  vertical  part  arises  from 

change  in  height  and  is  usually  smaller  than  the  horizontal  part. 

In  analogy  to  (21),  the  horizontal  part  is  given  by 
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£  denoting  the  dynamical  deflection  of  the  vertical  at  the  earth's 
surface;  cf.  (19)  and  Fig.  1.  For  the  vertical  part  we  have  by  (6); 

li  1_/I\  -1  /H  lllYr  (26) 

ah  3h  \y  )  ~  Y  \ah  "  y  3h/ 

or 

a <;  _  _  A£  _  (27 ) 

ah  r  y 

according  to  the  fundamental  equation  of  physical  geodesy  (PG,  p.  298,  eq. 

8-20)). 

Hence  (24)  becomes 


d<;  =  -eds  -  dh 


(28) 


On  integrating  this  relation  we  get  the  difference  of  the  height  anomaly 


The  gravity  anomaly  Ag  refers  to  the  earth's  surface  according  to  (8). 
The  first  term  on  the  right-hand  side  represents  the  Helmert  integral  (22) 
of  the  surface  deflection  £  ,  and  the  second  term  is  Mol odensky ' s 

correction  to  the  Helmert  integral,  necessary  to  obtain  height  anomalies. 
This  correction  depends  on  the  gravity  g  at  the  earth's  surface. 


4.  Reduction  for  Plumb-Line  Curvature 


For  Helmert's  formula  (22),  the  deflection  component  ec  refers  to 
the  geoid.  In  fact,  by  (13)  and  (14)  we  have 

(3 

eo  =  ^ o  ~^cosa  +  -^)cos4sina  , 

where  the  astronomical  coordinates  4  ,  A  are  taken  at  the  geoid  and 

o  o 

given  by  (10): 


4  =  $  +  6$ 

o 


A  =  A  +  &A  . 
o 


(31) 


Thus  the  astronomically  measured  surface  values  4  ,  A  must  be  reduced 

to  sea  level  by  applying  corrections  64  ,  6A  for  plumb-line  curvature. 

These  corrections  may  be  expressed  by 

,  iAcos*  =  t  (32 

where  OC  denotes  the  orthometric  correction  in  leveling  and  the  local 
axes  x  and  y  are  horizontal  and  directed  towards  north  and  east,  res¬ 
pectively.  Thence  follows  (PG,  p.  195): 


64 


.  H|a,!ta  tan6 
g  3x  ?  1 


5Acos4  =  -  3;  tanBo 

9  y  9 


(33) 


Here,  H  denotes  the  orthometric  height  (the  length  of  the  curved  plumb 
line  segment  PQP  in  Fig.  1),  g  is  the  mean  value  of  gravity  along  the 
plumb  line  between  PQ  and  P  ,  and  Si  and  Z2  designate  the  an¬ 
gles  of  inclination  of  a  terrain  profile  in  a  north-south  and  an  east-west 
direction. 
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Using  these  formulas,  we  thus  have  first  to  differentiate  the  ortho- 
metric  correction  in  a  horizontal  direction  (by  (32))  and  then  to  integrate 
by  the  Helmert  formula  (22).  As  a  matter  of  fact,  differentiation  and  sub¬ 
sequent  integration  should  give  back  the  original  quantity,  but  they  will 
in  general  fail  to  do  so  exactly  because  of  inaccuracies  inherent  in  the 
processes  of  numerical  differentiation  and  integration. 

Thus  it  is  preferable  to  apply  the  orthometric  correction  directly  to 
the  geoidal  height  difference:  by  PG,  pp.  200-201  we  have 


where  0C/\g  denotes  the  orthometric  correction  along  the  profile  AB 
and  is  the  "geometric  deflection  of  the  vertical"  at  the  earth's  surface 
given  by  (11)  and  (12);  cf.  a  corresponding  remark  at  the  end  of  Sec.  1. 

Eq.  (34)  thus  represents  a  classical  analogue  of  the  Molodensky 
formula  (29).  This  analogy  is  particularly  conspicuous  if  we  write  (34)  in 
differential  form,  using  (32)  and  (33): 

dN  =  -eds  +  ^  dg”  -  dH  . 

g  g 

The  comparison  of  this  expression  with  its  analogue  (28)  shows  as  the 
main  difference  the  fact  that  (35)  contains  mean  gravity  g"  .  Now, 
g  is  the  average  of  all  values  of  gravity  along  a  plumb-line  between 
sea  level  and  earth's  surface,  and  gravity  inside  the  earth  cannot  be 
measured,  nor  can  it  be  computed  rigorously  because  the  rock  density  p 
inside  the  earth  is  unknown.  Thus  g  cannot  be  determined  with 
complete  rigor. 

This  is  the  point  where  Molodensky's  criticism  of  the  classical  theo¬ 
ries  of  physical  geodesy  enters.  This  criticism  is  fully  justified  from  a 
conceptual  point  of  view  and  has  been  extremely  fruitful  for  the  de¬ 
velopment  of  modern  theoretical  geodesy.  On  the  other  hand,  from  a  practi¬ 
cal  point  of  view,  we  may  say  that,  with  a  reasonably  realistic 
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density  model,  mean  gravity  Tf  and  hence  orthometric  heights  H 
and  geoidal  heights  N  can  be  computed  quite  well  with  satisfactory 
accuracy;  cf.  the  simple  estimates  in  PG,  p.  169. 


5.  Approximate  Reduction  for  Plumb-Line  Curvature 

A  sufficiently  precise  determination  of  plumb-line  curvature  according 
to  (33)  or  of  the  orthometric  correction  in  (34)  is  rather  cumbersome, 
however. 

An  approximate  procedure  (Elmiger,  1969;  Gurtner,  1978)  uses  an 
analogy  to  classical  gravity  reduction,  applied  to  the  direction  of  the 
gravity  vector:  gravity  reduction  is  applied  to  the  magnitude  g  of  the 
gravity  vector  g. 


earth's  surface 


sea  level 
(geoid) 


FIGURE  5.  The  geometry  in  gravity  reduction 
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Consider  first  the  gravity  reduction  of  Poincare-Prey  (PG,  p.  165); 
cf.  Fig.  5.  It  consists  of  the  following  three  steps: 

1.  The  topographic  masses,  i.e.,  the  masses  above  the  geoid,  are  re¬ 
moved  computationally  by  subtracting  its  attraction  Ay  from  the  ob¬ 
served  gravity  value  g  at  the  surface  point  P  . 

2.  The  reduced  gravity  value  g  -  Ay  so  obtained  at  P  is  trans¬ 
ferred  to  point  PQ  at  sea  level  by  adding  the  free-air  reduction 

F  . 

3.  The  topographic  masses  are  restored  by  adding  its  attraction 
Ay  at  PQ  . 

The  result  of  these  three  steps, 

a  =q_A  +F+A0  (36) 

so  a  T  T 

gives  actual  gravity  gQ  on  the  geoid.  A  weak  point  of  the  procedure 
(apart  from  errors  in  Ay  and  A°  due  to  imperfect  knowledge  of 
density  p  )  is  the  computation  of  the  free-air  reduction1  F  by  the 
formula 


F  = 


(37) 


replacing  the  vertical  gravity  gradient  of  the  actual  gravity  field  (after 
removal  of  the  topographic  masses)  by  the  normal  gravity  gradient  3y/3h. 

The  astronomical  coordinates  $  and  A  (defining  the  direction  of 
the  gravity  vector  £  )  can  be  treated  in  complete  analogy.  Let  £y  he 
the  £  component  of  the  deflection  of  the  vertical  at  P  as  computed  from 
the  topographic  masses  only.  Let  further  be  £y°  the  corresponding 
topographic  deflection  of  the  vertical  at  Pq  .  In  complete  analogy  to 
the  three  steps  mentioned  above  we  have  now: 


1)  So  called  because  after  the  removal  of  the  topographic  masses,  the 
point  P  lies  "in  free  air". 


P 


9 


9 


9 


9 
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measured  at  P  . 

1.  removal  of 

topography,  at  P  . 

2.  "free-air  reduction"  P  -»•  PQ 

=  normal  plumb-line  curvature 

3.  restoration  of 

topography,  at  P0  . . . 

The  result  of  these  three  steps, 

*o  =  $  -  C-p  +  <$4>  +  Cy  , 


$ 

54> 


(38) 


thus  gives  the  astronomical  latitude  at  the  geoid.  This  formula  is 
completely  analogous  to  (36). 

The  comparison  of  (38)  with  (10)  yields,  for  the  actual  plumb-line 
curvature  6$  ,  the  approximate  expression 


6*  =  -  CT  +  Cj  +6« 


(39) 


where  6c?  denotes  the  normal  plumb-line  curvature  (16).  By  (10),  (11), 
and  (13)  we  thus  get  for  the  vertical  deflection  at  the  geoid: 

£q  =  -  <?  =  $-  <?  +  6<?  =  £  +  6$  ,  (40) 

with  an  analogous  equation  for  the  component  n 

Within  the  accuracy  of  this  procedure  we  thus  obtain  actual  vertical 
deflections  on  the  geoid,  in  the  same  way  as  the  Poi ncare-Prey  reduction 
gives  actual  gravity  on  the  geoid.  Hence  we  may  apply  to  £q  ,  nQ  the 
Helmert  formula  (22)  to  get  differences  of  the  geoidal  height  N 

If  we  have  a  reasonably  realistic  density  model  for  the  topographic 
masses,  we  can  compute  £7  ,  177  and  £7  ,  177  with  an  accuracy 

which  might  be  satisfactory  for  many  purposes.  The  weak  point  of  the 
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procedure,  in  the  same  way  as  for  gravity  reduction,  is  the  computation  of 
the  "free-air  reduction"  of  the  plumb  line  by  applying  the  normal 
plumb-line  curvature. 

Physically  this  implies  the  assumption  that,  after  the  removal  of  the 
topographic  masses,  the  earth's  crust  is  "regularized"  to  such  an  extent 
that  its  external  gravity  field  will  then  be  approximately  equal  to  the 
normal  ellipsoidal  gravity  field.  This  assumption  would  mean  that  the 
co-geoid  of  the  Bougner  reduction  (PG,  sec.  3-3)  coincides  with  the 
reference  ellipsoid.  In  fact,  however,  these  co-geoid  heights  are,  by  an 
order  of  magnitude,  larger  than  the  actual  geoid  heights  (Helmert,  1884, 
pp.  354-57;  PG,  sec.  3-6).  It  is  true  that  these  large  cogeoidal  heights 
have  a  smooth  and  regional  behavior  and  thus  affect  the  local  plumb-line 
curvature  to  a  lesser  degree,  but  the  assumption  of  a  normal  plumb-line 
curvature  is  nevertheless  rather  questionable. 

Topographic-isostatic  reduction.  A  considerably  bptter 
"regularization"  of  the  earth's  crust  may  be  expected  from  the  application 
of  a  topographic-isostatic  reduction.  The  Targe  regional  features  of  the 
co-geoid  are  essentially  reduced  in  this  way.  The  important  step  is  to 
take  some  isostatic  model,  and  it  is  not  so  essential  which  isostatic  model 
is  taken;  a  conventional  Ai ry-Heiskanen  model  with  T  =  30km  (PG,  Sec. 

3-5)  may  be  satisfactory  in  many  cases. 

Physically,  a  topographic-isostatic  reduction  means  that  the 
topographic  masses  (above  the  geoid)  are  computationally  removed  and  used 
to  fill  the  mass  deficiencies  which  exist  below  sea  level  according  to  the 
theory  of  isostasy.  If  the  isostatic  compensation  were  exact,  the  earth’s 
crust  would  become  completely  homogeneous  after  such  an  isostatic 
reduction,  so  that  the  irregularities  of  the  geoid  would  disappear  after 
reduction:  the  isostatic  co-geoid  would,  theoretically,  coincide  with  the 
reference  ellipsoid.  In  reality,  of  course,  no  isostatic  model  exactly 
corresponds  to  the  actual  geological  situation,  so  that,  rather  than 
coinciding  with  the  ellipsoid,  the  isostatic  co-geoid  will  deviate  less  and 
more  smoothly  from  the  ellipsoid  than  the  geoid  does.  Thus  the 
topographic-isostatic  reduction  achieves  a  considerable  smoothing  of 
geoidal  heights  and,  especially,  of  deflections  of  the  vertical;  cf.  Sec. 
10. 
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For  the  present  purpose  it  is  of  particular  interest  that  also  the 
curvature  of  the  plumb-line  will  be  smoothed  by  the  topographic-isostatic 
reduction.  In  this  way  it  will  certainly  become  closer  to  normal  plumb- 
line  curvature  than  with  a  mere  removal  of  the  topographic  masses. 

Therefore  it  is  appropriate  to  modify  the  reduction  procedure  given  at 
the  beginning  of  the  present  section,  in  such  a  way  that  in  Step  1  the 
topography  is  not  merely  removed  completely  but  that  it  is  transported  into 
the  interior  of  the  geoid  in  order  to  make  up  the  isostatic  mass  de¬ 
ficiencies.  This  procedure  has  to  be  reversed  in  Step  3  in  order  to  re¬ 
store  the  original  topography.  By  this  we  achieve  that  the  free-air  re¬ 
duction  by  means  of  the  normal  plumb-line  curvature  (Step  2)  appears  to  be 
better  justified. 

For  our  computation  this  means  that  the  vertical  deflections  Ct  » 

Oj  ,  computed  from  the  topographic  masses,  are  to  be  replaced  by  de¬ 
flections  Cjj  ,  n-fi  representing  the  combined  effect  of  the 
topographic  masses  and  their  isostatic  compensation.  The  numerical  com¬ 
putation  is  done  by  the  same  formulas  as  for  Sy  ,  dy  by  dividing 
the  topographic  and  compensating  masses  into  vertical  prisms,  if  possible 
using  a  digital  terrain  model1.  The  underlying  isostatic  model  should  be 
reasonably  realistic  and,  above  all,  simple  and  well  defined,  such  as  the 
Ai ry-Hei skanen  model. 

Then  eq.  (39)  is  to  be  replaced  by 

_  r  ,  rO  ,  -  ,  (41) 

<5$  =  ~  <iyj  +  £yj  +  6<J>  , 


1)  Relevant  formulas  can  be  found  in  (Forsberg  and  Tscherning,  1981) 
with  references  to  (MacMillan,  1958,  pp.  78-79)  and  (Jung,  1961, 
p.  167). 
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where  Sjl  and  5ti  express  the  deflection  effect  of  the  topographic 
and  compensating  masses  at  points  P  and  PQ  .  By  (40)  we  then  have 


€0  =  C  -  Cji  +  £Tj  +  <54>  , 


no  =  n  '  nTI  +  nTI 


(42) 


as  the  normal  curvature  component  iS  zero. 

Let  us  repeat  once  more  that  £  ,  n  denote  the  "geometrical"  sur¬ 

face  deflections  (11)  and  that  50  ,  no  denote  the  deflections  of  the 
vertical  on  the  geoid  (and  not  on  the  isostatic  co-geoid!).  Eq.  (22)  will 
then  give  differences  of  geoidal  heights  N  . 

It  should  be  noted,  however,  that  even  (41)  is  only  an  approximate  ex¬ 
pression  for  the  actual  plumb-line  curvature,  an  expression  which  is  better 
than  (39)  but  nevertheless  affected  by  a  certain  error. 

From  a  conceptual  point  of  view  let  us  emphasize  once  more  that  the 
geoidal  values  £0  ,  nQ  as  computed  from  (42)  correspond  to  actual 
gravity  anomalies  gQ  -  y  at  the  geoid,  where  g0  denotes  actual  gravity 
at  the  geoid  (inside  the  masses)  as  provided  by  the  reduction  of  Poincar^ 
and  Prey  (PG,  p.  146).  Topographic-isostatic  reduction  has  only  been  an 
auxiliary  device  to  compute  £0  ,  nQ  in  a  better  way.  In  the  next  sec¬ 

tion  we  shall,  however,  use  isostatic  reduction  in  a  way  that  is  analogous 
to  classical  geoid  determination  by  gravity  reduction. 
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6,  Determination  of  the  Co-geoid 

The  application  of  topographic-isostatic  reduction,  described  in  the 
preceding  section,  admits  of  an  interesting  alternative  variant.  In  the 
preceding  section,  the  isostatic  co-geoid  was  mentioned  to  provide  a  con¬ 
ceptual  background,  but  it  was  not  used  explicitly  since  ,  nQ  were 

applied  to  determine  the  geoid  directly.  The  alternative  approach  to  be 
described  now  determines  first  the  isostatic  co-geoid,  from  which  the  geoid 
is  obtained  by  taking  into  account  the  indirect  effect.  This  approach 
fully  corresponds  to  classical  gravimetric  geoid  determination  by  isostatic 
gravity  reduction  (PG,  Secs.  3-6  and  8-2). 

In  the  three-step  approach  described  in  the  preceding  section  we  will 
only  perform  Step  1  --  removal  of  topography  plus  isostatic  compensation  in 
their  effect  at  P  —  and  Step  2  --  free-air  reduction  from  P  to  PQ  by 
applying  the  normal  plumb-line  curvature  $<J>  .  The  third  step  -- 

restitution  of  topography  —  will  be  omitted.  The  result  is  thus 

5o  3  5  -  ?TI  +  6*  •  no  =  n  -  "TI  ;  (43) 


it  gives  boundary  values  since  after  the  removal  of  topographic  masses,  sea 

level  represents  a  boundary  of  the  solid  earth.  More  exactly,  this  bound- 

c  c 

ary  surface  is  the  isostatic  co-geoid.  The  deflections  £0  ,  nQ  as 

given  by  (*3)  are  the  precise  equivalent  of  isostatic  gravity  anomalies 

£gc=  £gi  (PG,  p.  140)  and  likewise  refer  to  the  co-geoid. 
o  c  c 

Then  it  is  possible  to  apply  the  Helmert  formula  (22)  to  £0  ,  nQ 

to  obtain  differences  of  co-geoid  heights  Nc  : 


(44) 
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with 


ec  =  £C cosa  +  nCsina 
o  o  o 


(45) 


From  the  co-geoid  obtained  in  this  way  we  get  the  actual  geoid  by  ap¬ 
plying  the  indirect  effect;  cf.  PG,  p.  141.  Thus,  the  geoidal  height  N 
follows  from 


N  =  NC  +  6N 


(46) 


For  the  indirect  effect  6N  ,  an  application  of  the  Rruns  formula  (6) 
gives 


6N 


(47) 


where  T°i  denotes  the  potential  of  the  topographic  masses  together  with  ,  ' 

their  isostatic  compensation,  taken  at  the  geoidal  point  P  .  This 

o  < 

potential  is  computed  similarly  as  the  effect  of  the  topographic-i sostati c  9. — - 

masses  on  the  gravity  anomaly  and  on  the  deflection  of  the  vertical. 

If  we  have  applied  this  procedure  correctly,  we  should  theoretically 
get  the  same  result  for  the  geoidal  height  N  as  if  we  had  used  the 

Helmert  integration  of  the  vertical  deflections  c,Q  .  given  by  (42).  *  - 

Thus  the  indirect  effect  <5N  must  satisfy  the  relation 


For  the  present  method,  this  result  is  of  little  use  since  cN  is  computed  ’• 

directly  from  (47);  perhaps  it  could  be  used  for  checking  purposes. 

9 _ 


An  advantage  of  the  present  procedure  is  the  fact  that  the  isostati- 
cally  reduced  vertical  deflections  ’  no  ’  much  smaller  and 

smoother,  can  be  interpolated  better  than  the  surface  values  £  ,  n  and 
also  better  than  the  geoidal  values  £q  ,  nQ  •  Being  only  an  es¬ 
sentially  equivalent  modification  of  the  procedure  of  Sec.  5,  it  shares, 
however,  its  main  drawback:  the  application  of  the  normal  plumb-line 
curvature  reduction  64>  is  problematical. 
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7.  Topographic-Isostatic  Reduction  from  a  Modern  Point  of  View 

For  the  reasons  mentioned  at  the  end  of  the  preceding  section,  it  is 
natural  to  try  and  find  a  way  which  makes  use  of  the  clear  advantages  of 
the  topographic-isostatic  reduction  but  avoids  the  problems  inherent  in  a 
free-air  reduction  from  P  to  PQ  . 

For  this  purpose  let  us  once  more  consider  eqs.  (43)  but  interpret 
them  differently.  We  shall  put 


5 


c 


=  K  - 


Cy  j  +  <$<!> 


(49) 


By  means  of  (18)  this  may  be  written 


5C  =  5  -  5- 


TI 


(50) 


The  interpretation  of  (50),  however,  is  clear,  simple,  and  rigorous: 
from  the  dynamic  deflections  of  the  vertical  at  P  ,  which  are  the  vpry 
quantities  T  and  n  ,  we  subtract  the  effect  of  the  topographic- 
isostatic  masses,  Sti  ancl  on  ,  likewise  at  P  .  The  vertical  de¬ 
flections  so  obtained,  and  tf  ,  thus  do  not  really  refer  to  the 

(co-)geoid;  in  reality,  they  refer  to  the  earth's  surface! 
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But  what,  then,  means  the  normal  plumb-line  curvature  5<J>  in  (49)? 
Hoes  it  not  mean  a  reduction  from  the  earth's  surface  to  sea  level?  No,  in 
eqs.  (18)  it  only  denotes  the  transformation  between  the  geometrical  and 
the  dynamical  deflection  of  the  vertical,  both  referred  to  the  point  P  of 
the  earth's  surface.  This  is  also  clear  from  Fig.  1,  which  illustrates  the 
formula 


e  =  e  +  <5 


(51) 


extending  (18)  to  an  arbitrary  azimuth,  6  being  defined  by  (20). 

This  interpretation  of  (49)  or  (50)  as  isostatical ly  reduced  de¬ 
flections  of  the  vertical  at  the  earth's  surface  is  exact ,  whereas  the 

interpretation  of  (43)  as  deflections  at  the  co-geoid  was  only  approximate. 

c  c 

Therefore  we  now  have  written  K  ,  n  instead  of  our  former  notation 
c  c 

t,o  ,  nQ  .  This  is  the  desired  rigorous  interpretation  of  our 
i sostatical ly  reduced  vertical  deflections. 

This  interpretation  exactly  corresponds  to  the  modern  view  of  gravity 
reduction  according  to  the  theory  of  Molodensky  as  described  in  PG,  sec. 
8-11.  According  to  this  view,  the  i sostati cal ly  (or  in  some  other  way) 
reduced  gravity  anomalies  continue  to  refer  to  the  earth's  surface.  The 
classical  gravity  reduction  (PG,  sec.  8-2)  had  comprised  two  procedures: 
mass  transport  and  shift  P  -*•  PQ  ;  the  new  view  of  gravity  reduction  only 
considers  the  mass  transport;  the  problematic  shift  P  -*•  PQ  is  avoided. 
Formally,  a  "normal  free-air  reduction": 

F=-|fh  (52) 

ah 


may  he  said  to  occur  also  in  Molodensky's  theory:  normal  gravity  Y  in 
the  new  definition  (8)  of  the  gravity  anomaly,  where  it  refers  to  the 
telluroid  point  0  ,  is  computed  by 
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Y  = 


+  lx 


3h 


h 


(53) 


with  h  =  QqQ  denoting  the  normal  height  of  P  .  But  instead  of  re¬ 
ducing  actual  gravity  g  downward ,  from  P  to  P  ,  now  normal  gravity 

is  reduced  upward.  Whereas  for  the  first  process  the  use  of  the  normal 
gradient  3y/ 3h  is  problematic,  it  is  fully  justified  for  the  second 
process. 

In  a  similar  way,  we  might,  if  we  wish,  interpret  6<j>  as  a  reduction 
of  4>  for  normal  curvature  of  the  plumb-line  upwards ,  say,  from  Pq  to 
P  .  This  is  possible  because  in  (15),  4>  cou  1  d  be  said  to  refer  to  P  ) 
(because  PQ  and  P^  practically  coincide),  and  because  <p  denotes  the 
latitude  of  the  tangent  to  the  normal  plumb  line  at  P  .  This  inter¬ 
pretation  is  instructive  because  of  the  analogy  with  gravity  reduction, 
though  regarding  <j>  and  as  geometric  and  dynamic  latitude  of  the  same 
point  P  appears  more  natural. 

As  pointed  out  above,  the  present  interpretation  of  £c  and  r,c  as 
isostatical ly  reduced  deflections  of  the  vertical  at  the  earth's  surface  is 
conceptually  rigorous  and  therefore  also  practically  more  accurate,  hut 
this  decisive  advantage  implies  a  computational  drawback  if  integration 
along  a  profile  is  used;  since  this  integration  must  now  be  performed  along 
the  earth's  surface  and  not  along  a  level  surface  such  as  the  geoid,  com¬ 
putation  will  be  more  complicated.  Instead  of  the  simple  Helmert  formula 
(22)  we  now  must  use  the  Molodensky  formula  (2°>): 


with 


c 

£  = 


c  c  . 

t,  cosa  +  n  sina 


(55) 
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gc  being  the  isostatical 1y  reduced  surface  value  of  gravity  (measured 

value  g  minus  attraction  of  the  topographic-isostatic  masses). 

c 

From  the  isostatic  height  anomalies  £  obtained  in  this  way,  we  then 
get  the  actual  height  anomalies  5  by  applying  the  indirect  effect: 


c  ^  , 

C  =  C  +  6; 


(56) 


with 


SC  =  hi  (57) 

Y 

This  is  completely  analogous  to  (46)  and  (47),  but  now  Tyi  is  the 
potential  of  the  topographi  .-i sostatic  masses  at  the  surface  point  P  . 

As  a  matter  of  fact,  normal  gravity  in  (47)  refers  to  the  ellipsoid,  and  in 
(57),  to  the  telluroid,  but  the  difference  is  generally  small. 

For  higher  mountains,  the  isostatic  reduction  procedure  described  in 
the  present  section  is  preferable  in  practice  to  a  direct  application  of 
Molodensky's  formula  (29)  because  the  i sostatical ly  reduced  vertical  de¬ 
flections  are  much  smoother  and  easier  to  interpolate.  It  is,  however,  ex¬ 
tremely  laborious  from  a  computational  point  of  view  since  the  integration 
must  be  performed  along  the  earth's  surface  (or,  what  is  practically  the 
same,  along  the  telluroid). 

The  procedure  described  in  Sec.  6  is  easier  and  less  laborious  since 

the  integration  is  performed  along  a  level  surface.  It  is,  however,  not 

rigorous  theoretically  and  less  accurate  practically.  Now  we  can  also 

understand  the  error  inherent  in  the  procedure  of  Sec.  6:  the  isostatic- 

c  c 

ally  reduced  deflections  ,  nQ  at  sea  level  are  by  (43)  simply  put 

equal  to  the  correspondi ng  surface  deflections  (49).  In  reality,  however, 

,  nC  and  are  related  by  analytical  continuation  and  by  no 
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means  identical.  For  the  concept  of  analytical  continuation  cf.(  PG,  sec. 
8-10  and  p.  324);  for  the  magnitude  of  the  effect  of  analytical 
continuation  on  C  see  Sec.  10. 

Finally,  we  remark  that  the  computational  drawback  of  the  present 
method,  the  Molodensky  integration  along  the  earth's  surface,  can  he 
completely  avoided  if  we  perform  our  computations  in  space:  instead  of 
integrating  along  a  surface,  we  perform  collocation  in  space.  This  modern 
procedure,  to  be  described  in  the  next  section,  permits  a  simple  and 
computationally  convenient  use  of  surface  deflections  and  also  their 
combination  with  gravimetric  and  other  data. 


32 


8.  Least-Squares  Collocation 


The  principle  of  collocation  is  very  simple.  The  anomalous  potential 
T  outside  the  earth  is  a  harmonic  function,  that  is,  it  satisfies 
Laplace's  differential  equation 


AT 


2  2 
a  T  .  3  T  . 

2  2 
3x£  3y 


(58) 


An  approximate  analytical  representation  of  the  external  potential  T  is 
obtained  by 


T(P )  =  f(P)  = 


(59) 


that  is,  by  a  linear  combination  f  of  suitable  base  functions  ^  , 

<>2  ,  ...  ,  <f>q  with  appropriate  coefficients  ^  .  All  these 

are  functions  of  the  space  point  P  under  consideration. 

As  T  is  harmonic  outside  the  earth's  surface,  it  is  natural  to 
choose  base  functions  4^  which  are  likewise  harmonic,  so  that 


=  0 


(60) 


in  correspondence  to  (58). 

There  are  many  simple  systems  of  functions  satisfying  the  harmonicity 
condition  (60),  and  thus  we  have  many  possibilities  for  a  suitable  choice 
of  base  functions  .  We  might,  for  instance,  choose  spherical  harmon¬ 

ics  or  potentials  of  suitably  distributed  point  masses,  depending  on 
whether  we  emphasize  global  or  local  applications. 
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The  coefficients  b|<  may  be  chosen  such  that  the  given  obser¬ 
vational  values  are  reproduced  exactly,  for  instance,  all  deflections  of 
the  vertical  in  a  given  area.  This  means  that  the  assumed  approximating 
function  f  in  (59)  gives  the  same  deflections  of  the  vertical  at  the  ob¬ 
servation  stations  as  the  actual  potential,  and  hence  may  well  be  con¬ 
sidered  a  suitable  approximation  for  T  . 

Let  us  now  try  and  put  these  ideas  into  a  mathematical  form. 

Interpolation.  Let  errorless  values  of  T  be  given  at  q  spatial 
points  Pi  ,  ?2  »  ...  ,  Pq  ;  these  points  may  lie  on  the 
earth's  surface  or  in  space  above  the  earth’s  surface.  Me  put 

T(Pi)  ■  fi  ,  1  -  1,  2,  ....  q  ,  (61) 


and  postulate  that  in  approximating  T(P)  by  f(P)  ,  the  observations 
(61)  will  be  reproduced  exactly.  The  condition  for  this  is 

£Vk<P1>-T(P1>*'l  •  <62> 

whence  the  system  of  linear  equations, 

J^1kbk  'fi  Wlth  Aik  =  VPi>  (63) 

or  in  matrix  notation. 


A  b  =  f  . 


(64) 


If  the  square  matrix  A_  is  regular,  then  the  coefficients  b^  are  uni¬ 
quely  determined  by 


» 


This  model  is  suitable,  for  instance,  for  a  determination  of  the  geoid 
by  satellite  altimetry,  since  this  method,  rather  directly,  yields  geoidal 
heights  Nj  and  hence,  hy  Bruns'  theorem  (6),  T(P-j)  =  Yi N-j 
For  the  astrogeodeti c  geoid  determination  we  must  generalize  this  model, 
which  leads  us  to 

Collocation.  Here  we  wish  to  reproduce,  by  means  of  the  approximation 
(59),  q  measured  values  which  again  are  assumed  to  be  errorless  (this  as¬ 
sumption  is  not  essential  and  will  be  dropped  later).  These  measured 
values  are  assumed  to  be  so-called  linear  functionals  L^T  ,  L^T  , 

...  ,  LqT  of  the  anomalous  potential  T  . 

In  fact,  deflections  of  the  vertical, 

131  I3j[  (66) 

^  ~  "  y  3x  ’  n  "  "  Y  3y  ’ 


hut  also  gravity  anomalies, 


^  ■  -  It  - f T 


(67) 


are  such  linear  functionals;  here,  xyz  denotes  a  local  coordinate  system 
in  which  the  z-axis  is  vertical  upwards  and  the  x  and  y  axes  are 
directed  towards  north  and  east,  and  R  =  6371km  is  a  mean  radius  of  the 
earth.  Eq .  (66)  is  a  consequence  of  (25),  with  3s  =  3x  or  3y  ;  normal 
gravity  Y  may  be  considered  constant  with  respect  to  horizontal  de¬ 
rivation.  Eq .  (67)  is  the  well-known  fundamental  equation  of  physical 
geodesy  in  spherical  approximation  (PG,  p.  298  and  88);  the  equations  (66) 
and  (67)  refer  to  the  earth's  surface. 


V 


By  saying  that  deflections  of  the  vertical  and  gravity  anomalies  are 
linear  functionals  of  T  ,  v/e  simply  indicate  the  fact  that  T  ,  n  , 
Ag  depend  on  T  by  the  expressions  (66)  and  (67)  which  clearly  are 
linear ;  they  are,  of  course,  the  linear  terms  of  a  Taylor  expansion, 
neglecting  quadratic  and  higher  terms.  In  the  above  notation  L-jT  ,  the 
symbol  Lj  denotes,  for  instance,  the  operation 


L  .12- 
i  y  3x 


applied  to  T  at  some  point, 
Putti ng 


4f  ■  4T  1  4 


and  substituting  (59)  we  get 


t,Bikbk  ’  *i  with  Bik '  4*k 


denotes  the  number  obtained  by  applying  the  operation  to 
the  base  function  ^  ;  the  coefficient  obtained  in  this  way 

does  not  depend  on  the  measured  values.  Eq .  (70)  is  a  linear  system  of  q 
equations  for  q  unknowns,  which  is  quite  similar  to  (63).  This  method  of 
fitting  an  analytical  approximating  function  to  a  number  of  given  linear 
functionals  is  called  collocation  and  is  frequently  used  in  numerical 
mathemati cs . 

It  is  clear  that  interpolation  is  a  simple  special  case  of  col¬ 
location,  in  which 


Lif  =  f(Pi) 
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is  the  "evaluation  functional",  giving  the  value  of  f  at  a  point 
P-j  .  Thus  we  see  that  in  both  interpolation  and  collocation,  the 
coefficients  bk  require  the  solution  of  a  linear  system  of  equations 
(which  in  general  will  not  be  symmetric). 

Least-squares  interpolation.  Let  us  consider  a  function 

K  =  K(P,Q)  , 


in  which  two  points  P  and  0  are  the  independent  variables.  Let  this 
function  K  be 

-  symmetric  with  respect  to  P  and  Q  , 

-  harmonic  with  respect  to  both  points,  everywhere  outside  a  certain 
sphere,  and 

-  positive-definite 

(the  positive  definitiveness  of  a  function  is  defined  similarly  as  in  the 
case  of  a  matrix).  Then  the  function  K(P,Q)  is  called  a  (harmonic) 
kernel  function;  cf,  (Moritz,  1980,  p.  205).  A  kernel  function  K(P,Q) 
may  serve  as  "building  material"  from  which  we  can  construct  base 
functions.  Taking  for  the  base  functions  the  form 


<0k(P)  =  K(P,Pk) 


where  P  denotes  the  variable  point  and  P|<  is  a  fixed  point  in  space, 
we  obtain  least-squares  interpolation. 

This  name  originates  from  the  statistical  interpretation  of  the  kernel 
function  as  a  covariance  function;  then  least-squares  interpolation  has 
some  minimum  properties  (least  variance,  similarly  as  in  least-squares 
adjustment).  This  interpretation  is  not  essential,  however;  one  may  also 
work  with  arbitrary  analytical  kernel  functions,  considering  the  procedure 
as  a  purely  analytical  mathematical  approximation  technique.  Normally  one 
tries  to  combine  both  aspects  in  a  reasonable  way. 


r 


Substituting  (73)  into  (63)  we  get 


Aik  =  K(Pi,FV  =  Cik  ; 


this  square  matrix  now  is  symmetric  (in  the  general  case,  is  not 

symmetric!)  and  positive  definite  (because  of  the  corresponding  properties 
of  the  function  K(P,Q)  ).  Then  the  coefficients  b|<  follow  from  (65) 
and  may  be  substituted  into  (59).  With  the  notation 

MP)  a  K(p.pu)  3  C„l  •  (75 


the  result  may  be  written  in  the  form 


f(p)  =  tcpl  cp2 


['ll  jil2  • 
u2 1  L22  • 


C  l'1  Tf  1 

clq  f1 

C2q  2 


|_Cql  Cq2  CqqJ  [fqJ 

known  from  least-squares  interpolation  of  gravity  (PG,  p.  268). 

Least-squares  collocation.  Here  we  again  derive  the  base  functions 
from  a  kernel  function  K(P,Q)  ,  but  in  a  way  slightly  different  from 
(73):  we  put 

*k(R)  3  lJk(P.Q)  ,  (77: 

where  L^Q  means  that  the  functional  Lk  is  applied  to  the  vari¬ 
able  Q  ;  the  result  no  longer  depends  on  0  (since  the  application  of  a 
functional  results  in  a  definite  number).  Thus  we  must  in  (7(1)  put 


Btk  ■  lW!«P.Q>  ■  cik  , 


(78) 


t.  .. 


which  gives  a  matrix  which  again  is  symmetric.  Solving  (70)  for  b|< 
and  substituting  into  (59)  gives  with 


Vp>  ■  *  cPk 


(79) 


the  formula 


f(P)  =  [Cp-j  Cpg  •••  Cpq] 


^12  •’*  £lq 
l22  L2q 

-1 

*1 

• 

C„9  •••  c 
q2  qq 

\ 

(80) 


This  is  formally  the  same  expression  as  (76),  but  with  f-j  replaced  by 
i -j  and  with  "covariances"  C-j|<  and  Cpi  defined  by  (78)  and 

(79) . 

In  the  statistical  interpretation,  f(P)  is  an  optimal  estimate  (in 
the  sense  of  least  variance)  for  the  anomalous  potential  T  and  hence  for 
the  height  anomaly  C  =  T/y  on  the  basis  of  arbitrary  measuring  data.  For 
geoid  determination  in  mountain  areas,  relevant  measuring  data  primarily 
are  §  ,  n  ,  and  £g  .  The  covariances  and  Cp-j  are  given 

by  known  analytical  expressions  (Moritz,  1980,  sec.  15).  A  general  com¬ 
puter  program  for  collocation  is  described  in  (Sunkel,  1980). 

Least-squares  collocation  may  easily  be  generalized  to  observational 
data  affected  by  random  errors;  systematic  effects  may  also  be  taken  into 
consideration.  In  addition  to  the  estimated  quantities  (  f  in  our  pre¬ 
sent  case)  we  may  also  compute  their  standard  error  by  a  formula  similar  to 

(80) .  A  comprehensive  presentation  of  a  least-squares  collocation  may  be 
found  in  (Moritz,  1980). 


It  is  well  known  that  the  direct  interpolation  of  free-air  gravity 
anomalies  (which  essentially  are  surface  gravity  anomalies  (8))  in  high 
mountains,  e.g.  by  least-squares  interpolation,  leads  to  relatively  poor 
results,  because  of  the  correlation  of  the  free-air  anomalies  with  ele¬ 
vation  (PG,  sec.  7-10).  This  correlation  with  elevation  constitutes  a  con 
siderable  trend,  which  must  be  removed  before  the  interpolation.  Bouguer 
anomalies  take  care  of  the  dependence  on  the  local  irregul arities  of  ele¬ 
vation;  isostatic  anomalies  are,  in  addition,  also  largely  independent  on 
the  regional  features  of  topography  (PG,  p.  285). 

In  exactly  the  same  way  we  must  remove  the  main  trend  of  the  vertical 
deflections  £  ,  n  and  the  gravity  anomalies  Ag  by  an  isostatic  re¬ 

duction,  before  applying  collocation.  Thus  isostatic  reduction,  pragmatic 
ally  regarded  as  trend  removal,  is  essential  for  the  practical  application 
of  least-squares  collocation  in  mountainous  regions  (Forsberg  and 
Tscherning,  1981). 

Physically  speaking,  we  transport  the  topographic  masses  to  the  inter 
ior  of  the  geoid  in  such  a  way  that  the  isostatic  mass  deficiencies  are 
filled.  The  observation  point  P  remains  in  its  position  on  the  earth's 
surface.  In  this  way,  not  only  the  harmonic  character  of  the  anomalous 
potential  T  outside  the  earth's  surface  is  preserved,  but,  in  addition, 
the  computational  removal  of  the  topographic  masses  above  sea  level  makes 
the  function  T  harmonic  down  to  sea  level.  Hence,  the  collocation 
formula  (80)  can  be  applied  also  at  sea  level,  giving  co-geoid  heights 
Nc  .  By  applying  the  inverse  reduction  (the  indirect  effect)  to  the 
computed  height  anomalies  and  co-geoid  heights  Nc  we  get  actual 
C  and  N  .  It  can  be  expected  that  errors  in  the  isostatic  model  used 
(e.g.,  an  Ai ry-Hei skanen  model)  will  largely  cancel  in  this  combined  pro¬ 
cedure  of  reduction  and  "anti -reduction." 

The  procedure  is  theoretically  optimal  and  well  suited  for  computer 
use.  The  integrability  conditions,  which  in  Helmert  integration  are 
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represented  by  the  closures  of  the  individual  triangles,  are  automatically 
taken  into  account.  The  fact  that  the  deflections  of  the  vertical  are 
given  only  in  a  certain  region  has  the  effect  that  the  geoid  can  only  be 
computed  in  that  region.  Since  even  by  collocation,  differences  in  geoidal 
height  between  two  neighboring  stations  A  and  B  depend  essentially  only 
on  the  deflections  in  those  two  stations,  the  lack  of  data  outside  the  re¬ 
gion  under  consideration  will  hardly  cause  a  noticeable  distortion.  Note, 
however,  that  the  addition  of  a  constant  to  all  geoidal  heights  N  will 
not  affect  the  deflections  of  the  vertical;  hence,  astrogeodetic  data  de¬ 
termine  the  geoidal  heights  only  up  to  an  additive  constant.  This  constant 
may  be  chosen  such  that  the  average  value  of  the  computed  N  is  zero,  and 
the  result  of  collocation  comes  near  to  this  case. 

To  get  immediately  nearly  geocentric  geoidal  heights  it  is  appropriate 
to  take  into  consideration  a  global  trend  which  mainly  affects  C  and  N 
itself,  by  subtracting  the  effect  of  a  suitable  global  gravity  field,  say 
the  gravity  earth  model  (given  as  a  spherical-harmonic  expansion  up  to  de¬ 
gree  180°  x  180°)  of  Rapp  (1981),  following  Sunkel  (1983).  This  will  be 
described  in  the  next  section;  in  the  present  section  we  shall  limit 
ourselves  to  the  isostatic  reduction. 

Computational  procedure.  It  consists  of  the  following  steps: 

1.  Transformation  of  the  astrogeodeti c  surface  deflections  £  ,  n 

from  the  local  datum  used  for  the  geocentric  Geodetic  Reference  System  1980 
by  the  well-known  differential  formulas  of  Vening  Meinesz  (PG,  sec.  5-9). 
This  is  necessary  since  collocation  requires  a  reference  system  which  is  as 
realistic  as  possible. 

2.  Application  of  the  normal  plumb  line  curvature  (16)  to  the 

"geometric"  surface  deflections  £  ,  n  gives  the  "dynamic"  surface  de¬ 
flections  £  ,  "n  by  (15). 

3.  Computation  of  the  gravity  anomalies  Ag  ,  also  referred  to  the 
earth's  surface  according  to  (8). 

4.  The  topographic-isostatic  reduction  of  "£  ,  "n  ,  Ag  by  (50) 

and  PG,  eq .  (3-94)  gives  values  £°  ,  nc  .  Agc  which,  of 

course,  continue  to  refer  to  the  surface  point  P  . 


[ 
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5.  The  application  of  collocation  to  ,  rf  ,  £g  gives 

height  anomalies  $  and  co-geoid  height  N  ,  by  simply  varying  the 
elevation  parameter  (  H  and  zero  ,  respectively)  in  the  collocation 
program. 

6.  By  applying  the  indirect  effect  (57)  and  (47)  we  get  actual 
height  anomalies  5  and  geoidal  heights  N  . 


9  . 
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10.  Geoid  Computation  for  Austria 

Sunkel  (1983)  has  used  least-squares  collocation  to  calculate  the 
geoid  for  the  main  part  of  Austria.  In  addition  to  the  isostatic  reduction 
according  to  Ai ry-Heiskanen  (T  =  30km)  ,  he  also  removes  a  global  trend 

by  means  of  a  earth  gravity  model,  represented  by  a  spherical -harmoni c 
expansion  up  to  a  certain  degree  N  .  In  particular,  he  used  the  model  of 
Rapp  (1981)  with  N  =  180 

After  removing  the  topographic-isostatic  trend  Tjj  and  this  global 
trend  (  EM  denotes  earth  model),  there  remains  a  residual  anomalous 

potential  <ST  ,  given  by 


S[  =  T 


+  T 


TI 


(81) 


N 

Since  the  earth  model  potential  is  represented  by  a  spherical- 

harmonic  expansion  up  to  degree  N  ,  it  may  be  appropriate  to  consider, 
for  the  isostatic  reduction,  only  the  effect  for  degrees  N  >  180°  , 

replacing  Ty  j  by 


(TTI)N>180°  =  TTI  “  TTI 


(82) 


where  represents  a  spherical-harmonic  expansion  for  Ty j  truncated 

at  degree  N  =  180  .  This  explains  eq.  (81). 

The  observations  ^  =  (  C  ,  n  ,  Ag  )  ,  which  represent 
linear  functionals  L-jT  .are  reduced  in  the  same  way,  obtaining 


*i  ‘  LiTTI  "  LiTEM  +  LiTTI  =  Li6T 


(83) 
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Adding  the  earth  model  reduction  to  the  computational  procedure  outlined  at 
the  end  of  the  preceding  section,  we  thus  have  the  following  flow  diagram: 


(LJ) 


observations  referred  to 
Geodetic  Reference  System  1930 


reduction 
TI,  EM 


-L  (T  +  TN  -  TM  ) 

1 '  TI  EM  Tl' 


(L  «T) 


!  collocation 


inverse  reduction 
TI,  EM 


*  T  .  T'  ■  T" 

TI  EM  TI 


N  =  (T/y)  ,  C  =  (T/y). 

O  il 
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Data .  The  topography  in  Austria  is  rather  varied,  with  elevations  up 
to  3800 m  .  The  density  of  astrogeodeti c  stations  was  If)  to  20  km  ;  the 
total  number  of  deflections  data  used  was  521  .  No  gravity  anomalies 

were  used  in  this  first  computation. 

The  topographic-i sostatic  reduction  of  the  deflections  of  the  vertical 
was  made  using  a  rather  crude  digital  terrain  model  consisting  of  mean  ele¬ 
vations  for  20"  x  20"  rectangles.  It  has  been  obtained  by  digitizing  a 
map  1:500000  .  The  standard  error  of  this  model  is  on  the  order  of 
+0.00m  .  Investigations  have  shown  that,  in  spite  of  its  poor  accuracy, 

the  model  is  reasonably  adequate  for  reduction  of  deflections  of  the  verti¬ 
cal  (it  is  totally  inadequate  for  gravity!).  In  fact,  the  reduction  error 
for  5  ,  n  is  approximately  proportional  to  terrain  inclination;  it  is 

thus  very  small  if  the  deflection  station  is  situated  in  an  area  of 
inclination  zero.  This  is  the  case  not  only  if  the  station  lies  in  a 
horizontal  plane,  but  also  if  it  lies  on  the  top  of  a  mountain,  as  most  de¬ 
flection  stations  do. 

For  the  collocation  computation,  the  covariance  function  of  Jordan  and 
Heller  (1978)  was  used  for  reasons  of  simplicity.  The  parameters  of  this 
function  were  determined  from  the  data. 

Results.  It  turned  out  that  almost  all  of  the  signal  (T,  N,  £) 
comes  from  the  topographic-isostatic  model  and  the  N  =  180  gravity  model 
used.  This  part,  TI  +  EM  ,  lies  between  41.5m  and  47.5m  .  The 
contribution  of  collocation  (y“1  T)  lies  between  jT).5m  ,  after 
removal  of  a  pronounced  trend  on  the  order  of  3m  . 

The  efficiency  of  topographic-isostatic  reduction  can  also  be  seen 
from  the  fact  that  topographi c-i sostat ic  has  reduced  the  variance  of  the 
deflections  of  the  vertical  in  Austria  (the  square  of  the  average  size  of  £ 

p  O 

and  n  )  from  30  arcsecc  to  5  arcsec*1  . 

Of  considerable  interest  is  the  effect  of  analytical  continuation  on 
the  isostatical ly  (+EM)  reduced  anomalous  potential  T  .  It  is  ex¬ 
pressed  hy  the  difference  y_1  T  at  the  earth's  surface  minus  Y"1  T 
at  sea  level.  This  difference  reaches  a  maximum  of  13  cm  in  the  Central 
Alps  and  is  otherwise  positive  and  negative. 
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Of  even  greater  interest  is  the  difference  between  height  anomalies 
(  =  y"*T  at  the  earth's  surface)  and  geoidal  heights  N  (  =  y-1T  at 
sea  level)  .  The  maximum  of  35cm  for  C  -  N  is  reached  at  the 
Grossglockner  mountain  (H  =  3800m)  .  The  results  are  in  good  agreement 
with  the  approximate  formula  (PG,  p.  327) 

C  -  N  =  -(981gal ) _1  AgRH 

where  AgR  is  the  Rouguer  anomaly  in  gal  and  H  is  the  elevation  in  the 
same  units  as  S  and  N  . 

A  comprehensive  information  on  the  geoid  in  Austria,  with  regard  to 
both  observations  and  computations,  can  he  found  in  (Lichtenegger,  et  a  1 . , 
1983),  of  which  the  main  results  have  also  been  presented  in  English  at  the 
XVIII  General  Assembly  of  IUGG/IAG  in  Hamburg,  August  1983  (Bretterbauer, 
1983;  Erker,  1983;  Grasegger,  et  al.,  1983;  Haitzmann,  et  al.,  1983). 
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